Introduction

Transition from vegetative growth to reproductive growth is an extremely vital process in life cycle of higher plants, and this is realized by flowering. The correct structure of a floral organ determines the success of fruit formation and affects seed yield greatly in crop plants.

The rapeseed (Brassica napus L.) is an important oil crop globally. Besides edible oil, it also provides feed rich in protein for animals and raw materials for biodiesel production. It is an allotetraploid (AACC, 2n=38) derived from the natural corssing between progenitors B. rapa (AA, 2n=20) and B. oleracea (CC, 2n=18) (Cai et al. 2012). The three most important yield-influencing factors in rapeseed are seed weight, seed number per silique and silique number per plant (Xu et al. 2014). The seed weight, usually indicated by 1000-seed weight (TSW), has been extensively studied and > 130 loci have been identified (Fu et al. 2015; Zhao et al. 2016; Sun et al. 2018). For instance, Geng et al. (2016) crossed a large-seed line with a small-seed line and then used the doubled haploid (DH) population to identify a hot region associated with the TSW trait on Chromosome A09, followed by screening out four genes highly related to seed weight based on annotations; Sun et al. (2018) detected 25 QTLs for TSW and 152 for seven seed-shape traits by a recombinant inbred line (RIL). Additionally, the multi-locular Brassica juncea, presenting higher locule number in siliuqe than normal, increased yield per plant, which was considered due to the production of more seeds in a silique (Xiao et al. 2013, 2018). Tetra-locular trait in B. rapa (subspecies trilocularis) was also reported: Yadava et al. (2014) fine-mapped the locus and revealed it was due to a CLAVATA3 homologous gene. However, reports about variation in silique number are relatively rare.

We previously reported a multi-silique trait in rapeseed (Chai et al. 2019). Unlike the multi-locular trait mentioned above, our multi-silique rapeseed line “zws-ms” showed three siliques on each carpopodium, which was resulted from its aberrant floral structure: three pistils, including two extra-outgrowth pistils, in one floral organ (Jiang et al. 1998; Chai et al. 2019). We located the underlying loci on chromosome A09 and C08, respectively; and the differentially expressed genes (DEGs) between the two near-isogenic lines (NILs) zws-217 (normal silique) and zws-ms at budding stage were preliminarily analyzed. However, subsequent developmental stages were not investigated. Thus far, Zhu et al. (2019a) performed RNA-seq in a natural mutant of rapeseed exhibiting abnormal differentiation of stems, finding that genes participating in the shoot apical meristem (SAM) activity maintenance, cytokinin biosynthesis, and signal transduction displayed greatly variation at transcriptional level. Lee et al. (2018) elucidated physiological and molecular mechanisms underlying seed development of the tetralocular silique in B. rapa by using the dynamic transcriptome profiling, which involved differentially expressed metabolic genes and changes in transcript abundance of genes in different developing seeds (21, 28, 35, and 42 days after flowering), as well as seed storage compounds. In addition, RNA-seq was also used for dynamic study in microRNAs (miRNAs) and their target genes grapevine (Vitis vinifera), by investigating their variations systematically and spatiotemporally (Wang et al. 2014).

In order to further investigate the molecular mechanism underlying this multi-silique trait in rapeseed, we dynamically analyzed transcriptome data in multi-silique zws-ms and its NILzws-217at both flowering stage and green podding stage in this study. DEGs, including up-regulated genes and down-regulated genes were identified. Particularly, DEGs sharing the same expression tendencies at the two stages were analyzed. Annotations based on Gene Ontology (GO) and the Kyoto Encyclopedia of Genes or Genomes (KEGG) were then further analyzed. These results provide a deepened understanding for this multi-silique morphology in rapeseed.

Materials and Methods

Plant material and growth conditions

The homozygous multi-silique rapeseed material “zws-ms” and its near isogenic line (NIL) “zws-217” were developed at the Crop Research Institute, Sichuan Academy of Agricultural Sciences. The only difference between the two lines was the silique(s) morphology (Fig. 1). They were grown in an experimental field in the Xindu (altitude of 472 m, 30°47′10″N 104°12′12″E), located in the Sichuan Basin, with an annual average temperature of 16.2°C in humid subtropical monsoon climate.

RNA isolation and library preparation

Three zws-ms plants (samples T07, T08, and T09) and three zws-217 plants (T10, T11, and T12) were selected respectively at flowering stage; while three individuals (T13, T14, and T15) from zws-ms line and three individuals (T16, T17, and T18) from zws-217 line were sampled, respectively, at green podding stage. The total RNA was isolated using an RNA Isolation Kit (Tiangen, Beijing, China). The quality control of RNA was then performed on NanoDrop 2000 (Thermo Fisher Scientific, Waltham, USA) by estimating concentration at OD260/OD230. The sequencing libraries were generated using an RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, USA) according to the given instructions.

Transcriptome sequencing

The RNA-seq was performed on the Illumina HiSeq X-ten platform. The initially produced raw reads were processed to generate clean reads by using in-house Perl scripts to remove adapter sequences, reads containing poly-N, or low-quality reads. Furthermore, the Q30, GC-content, and sequence duplicationlevels of the clean reads were calculated. Clean reads were then aligned to the reference genome of Brassica napus “Darmor-bzh” (www.genoscope.cns.fr/brassicanapus/data/) by using Tophat2 tools, in order to screen out the reads with a perfect match or one mismatch for consequent investigation.

Differentially expressed gene (DEG) analysis

DESeq R package (Wang et al. 2010) was used to detect DEGs. The P value was adjusted by controlling the false discovery rate (FDR) described by Benjamini and Hochberg (1995). Genes with an adjusted log₂FC of 2 and a P value < 0.01 were then identified as DEGs.

Annotation of genes

GO database was used to annotate the genes. The GO enrichment of the DEGs was calculated by the GOseq R packages (Young et al. 2010). The KEGG database (http://www.genome.jp/kegg/) was used to explore the high-level functions and utilities of the biological system (Kanehisa et al. 2007). KOBAS software (Mao et al. 2005) was used to test the statistical enrichment of the DEGs in various KEGG pathways. We also compared genes with their orthologs from model plant Arabidopsis, based on the TAIR database (www.arabidopsis.org/).

Results

Generation of the transcriptome sequencing

We performed the transcriptome sequencing (RNA-seq) and the six samples at flowering stage (T07, T08 and T09 from multi-silique line zws-ms and T10, T11 and T12 from NIL zws-217) totally generated 62.44 Gb of raw data, with average Q30 value of 90.40% and GC content of 47.61%; while the six samples at podding stage (T13, T14 and T15 from zws-ms and T16, T17 and T18 from zws-217) produced 56.58 Gb raw data in total, with average Q30 value of 90.13% and GC content of 47.84% (Table 1). For each sample generated approximately 34.9 M and 31.6 M clean reads at average in T07~T12 group (flower) and T13~T18 group (pod), respectively (Table 1). The proportion of total reads mapped to the reference genome in T07~T12 group ranged from 72.91 to 75.09%; while it ranged from 71.30 to 74.57% (Table 2).

Table 1: Summary of the transcriptome sequencing data

Sample	Number of clean reads	Number of clean bases	GC content (%)	% ≥ Q30
T07	39,872,033	11,902,667,750	47.38	90.28
T08	29,835,132	8,894,586,194	47.66	91.44
T09	37,493,875	11,172,228,764	47.49	90.14
T10	37,185,677	11,091,447,160	47.37	90.28
T11	34,064,197	10,166,051,118	47.84	90.38
T12	30,883,341	9,214,076,026	47.91	89.86
T13	39,613,493	11,825,231,742	47.93	89.78
T14	26,988,155	8,030,791,800	47.80	90.05
T15	39,013,568	11,636,718,296	47.83	89.48
T16	29,122,944	8,688,230,392	47.41	91.36
T17	27,011,899	8,047,537,522	48.05	89.69
T18	28,023,293	8,348,109,052	48.01	90.43

Note: T07, T08, and T09: stamens and pistils from three zws-ms individuals at the flowering stage; T10, T11, and T12: stamens and pistils of three zws-217 individuals at the flowering stage; T13, T14 and T15: pods from three zws-ms individuals at the green podding stage; T16, T17 and T18: pods of three zws-217 individuals at the green podding stage

Identification of differentially expressed genes (DEGs)

Fig. 1: The multi-silique trait in zws-ms compared with its near isogenic line, normal zws-217

Fig. 2: The number of differentially expressed genes (DEGs) in multi-silique line zws-ms compared with zws-217, at flowering and green podding stages

Genes with expression level fold change > 4 (log₂FC >2) and FDR <0.001 were identified as DEGs. In the stamens and pistils at the flowering stage, 123 DEGs were identified, among which 70 were up-regulated in zws-ms and 53 were down-regulated (Fig. 2, Table S1). In the green pods, we found 268 DEGs totally. Compared with zws-217, 219 of the genes were up-regulated while 49 were down-regulated in zws-ms (Fig. 2, Table S2). Further analysis found that 77 DEGs shared the same expression in the two stages. Precisely, 46 of these genes were up-regulated in multi-silique zws-ms at both flowering and podding stages, while 31 genes were down-regulated at both stages. The expression level of these genes, including 12 new genes, represented stably line specific.

Annotations for the DEGs

GO terms were generally divided into three categories: biological processes (BP), cellular components (CC), and molecular functions (MF). We analyzed the GO terms of all the DEGs between zws-ms and zws-217. At the flowering stage (Fig. 3a): the BP terms with the highest levels of enrichment included “cellular process (GO: 0009987)” and “metabolic process (GO: 0008152)”, involved in 58 and 52 DEGs, respectively; in the CC category, most enriched terms were “cell part (GO: 0044464)”, “cell (GO: 0005623)” and “organelle (GO: 0043226)”; in the MF category, “catalytic activity (GO: 0003824)” and “binding (GO: 0005488)” accounted the top enriched terms, involved in 46and 36 DEGs. While at the green podding stage (Fig. 3b), which provided more DEGs: the BP terms with the highest levels of enrichment included “cellular process (GO: 0009987)”, “metabolic process (GO: 0008152)” and “single-organism process (GO: 0044699)”, involved in 130, 137 and 126 DEGs, respectively; in the CC category, the three most enriched terms were “cell part (GO: 0044464)”, “cell (GO: 0005623)” and “organelle (GO: 0043226)”, the same as those at flowering stage; the top two enriched terms were also the same to those in flowers, but with larger amount in gene numbers, that was, 142 DEGs in “catalytic activity (GO: 0003824)” term and 86 DEGs in “binding (GO: 0005488)” term. It could be seen that at both flowering and podding stages, the DEGs represented the highly similar top enriched GO terms.

Table 2: Mapped reads from the transcriptome sequencing data

Sample	Total reads	Mapped reads	Unique mapped reads	Multiple mapped reads	Reads mapped to '+'	Reads mapped to '-'
T07	79,744,066	58,732,384	52,322,869	6,409,515	27,498,178	27,527,387
T08	59,670,264	44,805,427	40,097,058	4,708,369	21,066,421	21,076,839
T09	74,987,750	54,766,420	48,293,165	6,473,255	25,431,146	25,461,399
T10	74,371,354	54,308,878	46,502,178	7,806,700	24,557,783	24,584,252
T11	68,128,394	50,765,850	42,596,102	8,169,748	22,570,858	22,584,713
T12	61,766,682	45,032,876	37,624,886	7,407,990	19,986,073	20,007,944
T13	79,226,986	57,391,532	47,797,521	9,594,011	25,507,679	25,527,520
T14	53,976,310	39,199,615	32,621,936	6,577,679	17,378,575	17,386,457
T15	78,027,136	55,632,827	48,605,556	7,027,271	25,698,426	25,704,707
T16	58,245,888	43,040,388	39,108,157	3,932,231	20,596,179	20,600,427
T17	54,023,798	39,689,203	33,286,563	6,402,640	17,732,478	17,738,347
T18	56,046,586	41,796,426	35,683,464	6,112,962	18,976,344	18,991,906

Note: T07, T08, and T09: stamens and pistils from three zws-ms individuals at the flowering stage; T10, T11, and T12: stamens and pistils of three zws-217 individuals at the flowering stage; T13, T14, and T15: pods from three zws-ms individuals at the green podding stage; T16, T17, and T18: pods of three zws-217 individuals at the green podding stage

cKEGG pathway enrichment analysis revealed that at flowering stage (Fig. 4a), the two enriched pathways with most genes were “phagosome (ko04145)” and “oxidative phosphorylation (ko00190)”. “Phagosome (ko04145)” involved 6 DEGs: BnaA09g55930D, BnaC08g29080D, BnaC08g35500D, BnaC08g35720D, BnaCnng75420D and Cole_newGene_1984; “oxidative phosphorylation (ko00190)” enriched 5 DEGs: BnaC08g29080D, BnaC08g35720D, BnaC08g36360D, BnaCnng75420D and Cole_newGene_1984. But “vitamin B6 metabolism (ko00750)” got the highest value of enrichment factor, at 12.8, involving only one gene: BnaC08g37550D. At green podding stage (Fig. 4b), “starch and sucrose metabolism (ko00500)”, “pentose and glucuronate interconversions (ko00040)” and “carbon metabolism (ko01200)” had most genes enriched, six for each pathway respectively. The “flavonoid biosynthesis (ko00941)”, involved 4 genes BnaA03g60670D, BnaC04g18950D, BnaC07g37670D and BnaC09g17150D obtained the highest value of enrichment factor, at 12.2.

Fig. 3: Gene ontology (GO) terms of the DEGs at flowering and podding stages

a: GO terms at flowering stage; b: GO terms at green podding stage

Overall, 60 of the stably line-specific 77 genes were annotated by GO or KEGG (Table S3). When we aligned the above-mentioned 77 DEGs to model Arabidopsis, we found that 56 of these genes had corresponding orthologs in Arabidopsis, 9 genes of them did not get corresponding genes and 12 of them were new genes (Table 3). The chromosomes that harbored the most stably line-specific DEGs were chromosome C08 and A09, where we had identified the two associated regions according to our previous investigation (Chai et al. 2019). That was, 23 and 11 stably line-specific DEGs on them (Table 3), respectively.

Discussion

The present study performed RNA-seq at flowering and green podding stage, dynamically analyzed the transcriptional level of genes in zws-ms, a rapeseed material with multi-silique morphology. Similar trait (increased number of pistils) was reported in wheat (Triticum aestivum, Yang et al. 2015; Wei 2017; Guo et al. 2019; Zhu et al. 2019b), rye (Secale cereal, Malyshev et al. 2001), alfalfa (Medicago sativa, Nair et al. 2008) and so on. However, there has been little report about this trait in rapeseed. Multi-locular trait was reported in Brassica plants B. juncea (Xu

Table 3: Ortholog information of the stably line-specific DEGs

Gene in rapeseed		Ortholog in Arabidopsis
Gene ID	Regulated in zws-ms	Gene ID	Description
BnaA03g35870D	up	AT3G57550	guanylate kinase (AGK2)
BnaA07g04500D	up	AT2G04900	unknown protein
BnaA08g02930D	up	AT1G48920	nucleolin like 1 (NUC-L1)
BnaA09g06740D	up	AT5G64180	unknown protein
BnaA09g26320D	up	AT3G54620	basic leucine zipper 25 (BZIP25)
BnaA09g44370D	down	AT1G19000	Homeodomain-like superfamily protein
BnaA09g45300D	down	AT1G15000	serine carboxypeptidase-like 50 (scpl50)
BnaA09g45310D	up	AT1G14990	unknown protein
BnaA09g45320D	down	AT1G14980	chaperonin 10 (CPN10)
BnaA09g45610D	up	AT1G14330	Galactose oxidase/kelch repeat superfamily protein
BnaA09g46080D	down	AT1G13780	F-box/RNI-like/FBD-like domains-containing protein
BnaA09g46290D	up	AT5G03495	RNA-binding (RRM/RBD/RNP motifs) family protein
BnaA09g48320D	down	AT1G09690	Translation protein SH3-like family protein
BnaA09g49480D	down	AT1G07390	receptor like protein 1 (RLP1)
BnaAnng30260D	up	AT3G54620	basic leucine zipper 25 (BZIP25)
BnaC01g02500D	up	AT5G35640	Putative endonuclease or glycosyl hydrolase
BnaC01g43270D	up	AT5G57590	biotin auxotroph 1 (BIO1)
BnaC02g06570D	down	AT5G16400	thioredoxin F2 (TRXF2)
BnaC03g23820D	up	AT2G42670	Protein of unknown function (DUF1637)
BnaC03g57080D	up	AT3G47630	CONTAINS InterPro DOMAIN/s: Mitochondrial matrix Mmp37
BnaC04g29730D	up	AT1G76680	12-oxophytodienoate reductase 1 (OPR1)
BnaC04g39120D	down	AT2G27450	nitrilase-like protein 1 (NLP1)
BnaC05g26860D	down	AT1G14800	Nucleic acid-binding, OB-fold-like protein
BnaC06g16950D	up	AT3G59000	F-box/RNI-like superfamily protein
BnaC07g33980D	up	AT4G16900	Disease resistance protein (TIR-NBS-LRR class) family
BnaC08g29060D	down	AT1G12820	auxin signaling F-box 3 (AFB3)
BnaC08g35720D	up	AT2G21410	vacuolar proton ATPase A2 (VHA-A2)
BnaC08g35850D	up	AT2G21300	ATP binding microtubule motor family protein
BnaC08g36200D	down	AT2G20820	unknown protein
BnaC08g37460D	down	AT1G17980	poly(A) polymerase 1 (PAPS1)
BnaC08g38300D	down	AT1G79100	arginine/serine-rich protein-related
BnaC08g39020D	up	AT1G15130	Endosomal targeting BRO1-like domain-containing protein
BnaC08g39120D	up	AT1G14990	unknown protein
BnaC08g39130D	up	AT1G14980	chaperonin 10 (CPN10)
BnaC08g39360D	up	AT1G14720	xyloglucan endotransglucosylase/hydrolase 28 (XTH28)
BnaC08g39990D	up	AT1G14000	VH1-interacting kinase (VIK)
BnaC08g40040D	down	AT1G13900	Purple acid phosphatases superfamily protein
BnaC08g40320D	up	AT1G13450	Homeodomain-like superfamily protein
BnaC08g40410D	up	AT5G19320	RAN GTPase activating protein 2 (RANGAP2)
BnaC08g40810D	up	AT1G12970	plant intracellular ras group-related LRR 3 (PIRL3)
BnaC08g41390D	down	AT1G12240	ATBETAFRUCT4
BnaC08g41540D	down	AT1G12140	flavin-monooxygenase glucosinolate S-oxygenase 5 (FMO GS-OX5)
BnaC08g41720D	down	AT1G11910	aspartic proteinase A1 (APA1)
BnaC08g42280D	down	AT1G10930	RECQ4A
BnaC08g42450D	down	AT1G09970	LRR XI-23
BnaC08g49500D	up	AT1G18720	Protein of unknown function (DUF962)
BnaC08g49610D	down	AT1G10760	STARCH EXCESS 1 (SEX1)
BnaC08g49940D	down	AT1G10720	BSD domain-containing protein
BnaC09g05960D	up	AT5G63550	DEK domain-containing chromatin associated protein
BnaC09g06220D	up	AT5G64090	FUNCTIONS IN: molecular_function unknown
BnaC09g06260D	up	AT5G64180	unknown protein
BnaCnng17490D	up	AT3G47680	DNA binding
BnaCnng23190D	down	AT1G16880	uridylyltransferase-related
BnaCnng24040D	up	AT1G53310	phosphoenolpyruvate carboxylase 1 (PPC1)
BnaCnng68410D	up	AT3G57290	eukaryotic translation initiation factor 3E (EIF3E)
BnaCnng75420D	up	AT2G21410	vacuolar proton ATPase A2 (VHA-A2)
BnaA04g06410D	up	-	-
BnaAnng13790D	up	-	-
BnaC03g09190D	down	-	-
BnaC03g18930D	up	-	-
BnaC03g19830D	up	-	-
BnaC04g10370D	down	-	-
BnaC06g07110D	up	-	-
BnaC06g42000D	up	-	-
BnaC07g36960D	up	-	-
Cole_newGene_1983	down	-	-
Cole_newGene_1984	down	-	-
Cole_newGene_2071	down	-	-
Cole_newGene_2073	up	-	-
Cole_newGene_2243	down	-	-
Cole_newGene_269	up	-	-
Cole_newGene_3682	down	-	-
Cole_newGene_4151	down	-	-
Cole_newGene_4384	up	-	-
Cole_newGene_4855	up	-	-
Cole_newGene_6687	up	-	-
Cole_newGene_8277	down	-	-

et al. 2014; Xiao et al. 2018) and B. rapa (Fan et al. 2014; Yadava et al. 2014; Lee et al. 2018). But that trait increased the seed number per silique, rather than silique number per plant, and tqas different from our zws-ms line.

Fig. 4: Statistics of KEGG Pathway Enrichment

a: KEGG enrichment at flowering stage; b: KEGG enrichment at green podding stage

This multi-silique trait was originally discovered in hybridization progeny of B. napus × B. rapa and in previous researches (Jiang et al. 1998; Chai et al. 2019), we performed the association analysis combined with bulked segregant analysis (BSA) and detected two associated regions on chromosome A09 and C08, respectively. After analyzing the transcriptional levels at budding stage and annotations of the genes, we screened out some potential candidate genes. In order to obtain a more comprehensive understanding of the variations in zws-ms, which might cause the difference to its NIL zws-27, we further performed RNA-seq at flowering and green podding stages, dynamically displaying the transcriptional level.

Transcriptome sequencing (RNA-seq) can rapidly detect all transcripts from a specific tissue or in a particular state, and it provides a general profile for transcriptional levels of all expressed genes, the survey of alternative splicing events and even synteny analysis between sub-genomes (Paritosh et al. 2014; Vitulo et al. 2014; Xie et al. 2015). In this study, RNA-seq generated approximately 9.9 Gb for each sample (including both flowering and green podding stage), this adequate data was then validated by the Q30 value and the proportion of reads mapped to the reference genome.

Difference in expression level was indicated by log₂FC and DEGs were identified. At flowering stage, 70 DEGs were found up-regulated and 53 DEGs were down-regulated in multi-silique zws-ms; at the podding stage, the number of total DEGs increased significantly, especially the up-regulated genes in zws-ms. GO enrichment analysis showed that data from the two stages were highly similar in the most enriched terms in each category. Some genes were annotated to terms concerning flower/fruit structure or development: BnaC08g29060D, showing no expression in zws-ms, was annotated to “photoperiodism, flowering (GO: 0048573)”; BnaC08g39360D was annotated to “fruit development (GO:0010154)”; the term “ovule development (GO:0048481)” involved three genes: BnaC08g41780D, BnaA09g45320D and BnaC08g39130D. Moreover, these genes were all on chromosome A09 or C08, where we previously identified two associated regions (Chai et al. 2019). As to the KEGG pathways, the enriched pathways at two stages were different, but it also found that in the most enriched pathways, the genes accounting highest proportion were also from chromosome A09 or C08. This highlights the potential importance of these genes.

Arabidopsis and B. napus have a common ancestor (Cai et al. 2012), and Arabidopsis is a well-studied model plants, of which most of the genes were investigated adequately. So orthologs from Arabidopsis provide a reference for understanding the functions of their corresponding rapeseed genes. Fifty six out of the 77 stably line-specific genes found orthologs in Arabidopsis. After we further analyzed them and combined with our earlier report (Chai et al. 2019), we screened out following genes: BnaC08g29060D, showing no expression in zws-ms, was annotated to “photoperiodism, flowering (GO: 0048573)”. The ortholog in Arabidopsis is AT1G12820, which is known as auxin signaling F-box 3 (AFB3). It is induced by nitrate in primary roots and the auxin receptor involved in primary and lateral root growth inhibition (Dong et al. 2006). AT1G14720 encodes the xyloglucan endotransglucosylase/hydrolase 28 (XTH28), a member of Glycoside Hydrolase Family 16. In Arabidopsis atxth28 mutant displayed abnormal floral structure: the lower elongation ratio of the stamen at the pollination stage, the increased angle between the stamen and pistil, and the aberrant orientation of anther dehiscence (Kurasawa et al. 2008). Thus, its ortholog BnaC08g39360D in rapeseed, is supposed to cause some unknown changes in flower potentially. The term “ovule development (GO:0048481)” involved three genes: BnaC08g41780D, BnaA09g45320D and BnaC08g39130D. In fact, these three genes have been highlighted in our earlier report (Chai et al. 2019): AT1G11870, homologous to BnaC08g41780D, is identified as OVULE ABORTION 7 (OVA7), of which the disruption could result in ovule abortion in Arabidopsis (Berg et al. 2005); thus the different expression of BnaC08g41780D between zws-ms and zws-217 is also considered to result in some unknown process in fruit development in rapeseed herein. In addition to “ovule development”, BnaA09g45320D and BnaC08g39130D were also annotated to “response to cold (GO:0009409)”. It is notable that although they shared the same ortholog (AT1G14980), the two genes represented completely opposite expression model. That is, BnaA09g45320D was downregulated while BnaC08g39130D was unregulated in zws-ms in both flower and pod, which is consistent with our earlier study in bud (Chai et al. 2019). Although how these hypothetic chaperonins, as well as those genes mentioned above, causes abnormal flower development needs further investigation, this study highlights their potential functions in this multi-silique trait, especially BnaC08g39360D, BnaC08g41780D, BnaA09g45320D and BnaC08g39130D.

Conclusion

We performed the RNA-seq to analyze the DEGs between multi-silique rapeseed zws-ms and its NIL zws-217. 123 and 268 DEGs were identified at the flowering and green podding stage, respectively. Further analysis selected out 77 line-specific genes, which were expressed only in either zws-ms or zws-217. Combined with the information conferred by GO annotations, KEGG pathways and ortholog information from Arabidopsis, we identified some potential candidate genes underlying the multi-silique trait, like BnaC08g39360D, BnaC08g41780D, BnaA09g45320D and BnaC08g39130D. This deepens the understanding of molecular mechanism for the multi-silique phenomenon in rapeseed.

Acknowledgments

We thank the following funding for supporting this study: the National Key Research and Development Plan- International Cooperation Plan of Ministry of Science and Technology (2018YFE0108200); the International Cooperation Plan of SAAS (CGZH2019GH01); the Modern Agro-Industry Technology Research System of China (CARS-12); Young Leading Talents Project of Sichuan Academy of Agricultural Sciences (2019LJRC003 and 2019LJRC002).

Author Contributions

LC and JZ planned and carried out the experiments, performed data analysis, and drafted the manuscript, HL, JJ, CC, and BZ revised the manuscript and optimized the figures; PH and LZ assisted with DNA extraction and data analysis; LJ designed and supervised the experiments.

References

Benjamini Y, Y Hochberg (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. J Roy Stat Soc B 57:289‒300

Berg M, R Rogers, R Muralla, D Meinke (2005). Requirement of aminoacyl-tRNA synthetases for gametogenesis and embryo development in Arabidopsis. Plant J 44:866‒878

Cai G, Q Yang, Q Yang, Z Zhao, H Chen, J Wu, C Fan, Y Zhou (2012). Identification of candidate genes of QTLs for seed weight in Brassica napus through comparative mapping among Arabidopsis and Brassica species. BMC Genet 13:105-121

Chai L, J Zhang, K Lu, H Li, L Wu, H Wan, B Zheng, C Cui, J Jiang, L Jiang (2019). Identification of genomic regions associated with multi-silique trait in Brassica napus. BMC Genomics 20:304-322

Dong L, L Wang, YE Zhang, Y Zhang, X Deng, Y Xue (2006). An auxin-inducible F-box protein CEGENDUO negatively regulates auxin-mediated lateral root formation in Arabidopsis. Plant Mol Biol 60:599‒615

Fan C, Y Wu, Q Yang, Y Yang, Q Meng, K Zhang, J Li, J Wang, Y Zhou (2014). A novel single-nucleotide mutation in a CLAVATA3 gene homolog controls a multilocular silique trait in Brassica rapa L. Mol Plant 7:1788‒1792

Fu Y, D Wei, H Dong, Y He, Y Cui, J Mei, H Wan, J Li, R Snowdon, W Friedt, X Li, W Qian (2015). Comparative quantitative trait loci for silique length and seed weight in Brassica napus. Sci Rep 5; Article 14407

Geng X, C Jiang, J Yang, L Wang, X Wu, W Wei (2016). Rapid identification of candidate genes for seed weight using the SLAF-Seq method in Brassica napus. PLoS One 11; Article e0147580

Guo J, G Zhang, Y Song, Z Li, S Ma, N Niu, J Wang (2019). Comparative proteomic analysis of multi-ovary wheat under heterogeneous cytoplasm suppression. BMC Plant Biol 19; Article 175

Jiang LC, QX Zhang, XB Pu, ZY Zhang (1998) Evaluation of genetic and physiological Characteristics of aggregate-siliqua rapeseed. Southwest Chin J Agric Sci 11:62‒68 (in Chinese)

Kanehisa M, M Araki, S Goto, M Hattori, M Hirakawa, M Itoh, T Katayama, S Kawashima, S Okuda, T Tokimatsu, Y Yamanishi (2007). KEGG for linking genomes to life and the environment. Nucl Acids Res 36:480‒484

Kurasawa K, A Matsui, R Yokoyama, T Kuriyama, T Yoshizumi, M Matsui, K Suwabe, M Watanabe, K Nishitani (2008). The AtXTH28 Gene, a xyloglucan endotransglucosylase/hydrolase, is involved in automatic self-pollination in Arabidopsis thaliana. Plant Cell Physiol 50:413‒422

Lee Y, K Kim, J Lee, Y Cha, Y Moon, Y Song, E Jeong, S Ahn, W Park (2018). Comprehensive transcriptome profiling in relation to seed storage compounds in tetralocular Brassica rapa. J Plant Growth Regul 37:867‒882

Malyshev S, V Korzun, A Voylokov, V Smirnov, A Börner (2001). Linkage mapping of mutant loci in rye (Secale cereale L.). Theor Appl Genet 103:70‒74

Mao X, T Cai, JG Olyarchuk, L Wei (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21:3787‒3793

Nair RM, DM Peck, IS Dundas, DA Samac, A Moore, JW Randles (2008). Morphological characterisation and genetic analysis of a bi-pistil mutant (bip) in Medicago truncatula Gaertn. Sex Plant Reprod 21:133‒141

Paritosh K, V Gupta, SK Yadava, P Singh, AK Pradhan, D Pental (2014). RNA-seq based SNPs for mapping in Brassica juncea (AABB): Synteny analysis between the two constituent genomes A (from B. rapa) and B (from B. nigra) shows highly divergent gene block arrangement and unique block fragmentation patterns. BMC Genomics 15; Article 396

Sun L, X Wang, K Yu, W Li, Q Peng, F Chen, W Zhang, S Fu, D Xiong, P Chu, R Guan, J Zhang (2018). Mapping of QTLs controlling seed weight and seed-shape traits in Brassica napus L. using a high-density SNP map. Euphytica 214:228–240

Vitulo N, C Forcato, E Carpinelli, A Telatin, D Campagna, M D'Angelo, R Zimbello, M Corso, A Vannozzi, C Bonghi, M Lucchin, G Valle (2014). A deep survey of alternative splicing in grape reveals changes in the splicing machinery related to tissue, stress condition and genotype. BMC Plant Biol 14; Article 99

Wang C, X Leng, Y Zhang, E Kayesh, Y Zhang, X Sun, J Fang (2014). Transcriptome-wide analysis of dynamic variations in regulation modes of grapevine microRNAs on their target genes during grapevine development. Plant Mol Biol 84:269‒285

Wang L, Z Feng, X Wang, X Wang, X Zhang (2010). DEGseq: An R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26:136‒138

Wei SH (2017). Characterization and expression of WAG-2 transcripts in a wheat three-pistil mutant line. Russ J Plant Physiol 64:680‒687

Xiao L, X Li, F Liu, Z Zhao, L Xu, C Chen, Y Wang, G Shang, D Du (2018). Mutations in the CDS and promoter of BjuA07.CLV1 cause a multilocular trait in Brassica juncea. Sci Rep 8:5339-5351

Xiao L, H Zhao, Z Zhao, D Du, L Xu, Y Yao, Z Zhao, X Xing, G Shang, H Zhao (2013). Genetic and physical fine mapping of a multilocular gene Bjln1 in Brassica juncea to a 208-kb region. Mol Breed 32:373‒383

Xie BB, D Li, WL Shi, QL Qin, XW Wang, JC Rong, CY Sun, F Huang, XY Zhang, XW Dong, XL Chen, BC Zhou, YZ Zhang, XY Song (2015). Deep RNA sequencing reveals a high frequency of alternative splicing events in the fungus Trichoderma longibrachiatum. BMC Genomics 16; Article 54

Xu P, Z Lv, X Zhang, X Wang, Y Pu, H Wang, B Yi, J Wen, C Ma, J Tu, T Fu, J Shen (2014). Identification of molecular markers linked to trilocular gene (mc1) in Brassica juncea L. Mol Breed 33:425‒434

Yadava SK, K Paritosh, P Panjabi-Massand, V Gupta, A Chandra, YS Sodhi, AK Pradhan, D Pental (2014). Tetralocular ovary and high silique width in yellow sarson lines of Brassica rapa (subspecies trilocularis) are due to a mutation in Bra034340 gene, a homologue of CLAVATA3 in Arabidopsis. Theor Appl Genet 127:2359‒2369

Yang Z, Z Peng, S Wei, M Liao, Y Yu, Z Jang (2015). Pistillody mutant reveals key insights into stamen and pistil development in wheat (Triticum aestivum L.). BMC Genomics 16; Article 211

Young MD, MJ Wakefield, GK Smyth, A Oshlack (2010). Gene ontology analysis for RNA-seq: Accounting for selection bias. Genome Biol 11; Article R14

Zhao W, X Wang, H Wang, J Tian, B Li, L Chen, H Chao, Y Long, J Xiang, J Gan, W Liang, M Li (2016). Genome-wide identification of QTL for seed yield and yield-related traits and construction of a high-density consensus map for QTL comparison in Brassica napus. Front Plant Sci 7; Article 17

Zhu KM, S Xu, KX Li, S Chen, S Zafar, W Cao, Z Wang, LN Ding, YH Yang, YM Li, XL Tan (2019a). Transcriptome analysis of the irregular shape of shoot apical meristem in dt (dou tou) mutant of Brassica napus L. Mol Breed 39:39–49

Zhu X, Y Ni, R He, Y Jiang, Q Li, J Niu (2019b). Genetic mapping and expressivity of a wheat multi-pistil gene in mutant 12TP. J Integr Agric 18:532‒538